%--------------------------------------------------------------------------
function [x0, tf] = continuation_dro_(x0,auxdata)
mu = auxdata.mu;
RelTol = auxdata.RelTol;
AbsTol = auxdata.AbsTol;

MaxNum = 100; % 最大迭代次数
Maxtf = 200; % 最长积分时间
DAMP = 0.98; % 增加DAMP改善收敛性

% -------------------------------------------------------------------------
xf = ones(6,1);
count = 1;		 % 开始计数

while (abs(xf(1) - x0(1)) > 1e-14)
    % 超出最大迭代次数
    if count > MaxNum
        disp('Maximum iterations exceeded') ;
        break
    end
    
    % 状态变量和STM
    xx0 = zeros(42,1);
    xx0(1:6) = x0;
    xx0(7:42) = reshape(eye(6) , 36 , 1);

    %求解ode
    options = odeset('RelTol',RelTol,'AbsTol',AbsTol,'Events',@event_);
    [tt,xx,tf,xxf,ie] = ode113(@fun_crtbp_stm_,[0,Maxtf],xx0,options,auxdata);
    xf = xxf(2,1:6); % 终止状态变量
    phif = reshape(xxf(2,7 : end),6,6); % 单值矩阵

    %更新相对距离
    d2 = (xf(1) + mu)^2 + xf(2)^2 + xf(3)^2;
    r2 = (xf(1) - (1 - mu))^2 + xf(2)^2 + xf(3)^2;
    d3 = d2^1.5;
    r3 = r2^1.5;
    
    %更新势能
    Ux = xf(1)  - (1 - mu) * (xf(1) + mu) / d3 - mu * (xf(1) - (1 - mu)) / r3;
    Uy = xf(2) - (1 - mu) * xf(2) / d3 - mu*xf(2) / r3;
    
    %dxf和yf的二阶导
    ddxf = Ux + 2*xf(5);
    ddyf = Uy - 2*xf(4);
    
    %求解初始修正量
    delta_x0 = pinv(phif(1,5))*(xf(1)-x0(1));
    %delta_x0 = pinv([phif(1,5);phif(4,5);phif(5,5)-1]'*[phif(1,5);phif(4,5);phif(5,5)-1])*[phif(1,5);phif(4,5);phif(5,5)-1]'*[xf(1)-x0(1) ; xf(4) ; xf(5)-x0(5)];
    x0(5) = x0(5) - (1 - DAMP^count) * delta_x0(1);
   
    
    count = count + 1;
    %DAMP = 0.99 / sqrt(count); % 随着迭代次数增加，逐渐减小阻尼因子
    
   
end

end
% -------------------------------------------------------------------------

% 积分终止条件，即积分到x-z平面终止
function [value,isterminal,direction] = event_(t, xx, auxdata)
    value = xx(2); % 检测 x(2) = 0 
    isterminal = 1; % 终止积分
    direction = 1; % 正向穿过
end